Analytic, Group-Theoretic Density Profiles for Confined, Correlated A^-Body Systems 
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Confined quantum systems involving A'^ identical interacting particles are to be found in many 
areas of physics, including condensed matter, atomic and chemical physics. A beyond-mean-field per- 
turbation method that is applicable, in principle, to weakly, intermediate, and strongly-interacting 
systems has been set forth by the authors in a previous series of papers. Dimensional perturbation 
theory was used, and in conjunction with group theory, an analytic beyond-mean-field correlated 
wave function at lowest order for a system under spherical confinement with a general two-body 
interaction was derived. In the present paper, we use this analytic wave function to derive the 
corresponding lowest-order, analytic density profile and apply it to the example of a Bose-Einstein 
condensate. 

PACS numbers: 03.65.Ge,31.15.Hz,31.15.Md,03.75.Hh 



I. INTRODUCTION 

Confined quantum systems are widespread across 
many areas of physics. They include, among 
other examples, quantum dots[l], two-dimensional elec- 
tronic systems in a corbino disk geometry [2], Bose- 
Einstein condensates (BECs)[3], rotating superfluid he- 
lium systems [4], and atoms (where the massive nucleus 
provides the confining potential for the electrons). These 
new (and sometimes old) systems possess from a few tens 
to millions of particles and present a challenge for ex- 
isting A^-body methods when the mean-field approach 
fails [5]. 

These systems span an enormous range in interparticle 
interaction strength. For example, a gaseous BEG is typ- 
ically a weakly-interacting system for which a mean-field 
description is perfectly adequate. On the other hand, 
a superfluid helium system is a strongly interacting sys- 
tem. In this regard we note that the scattering length 
of a gaseous BEG tuned by a magnetic field can cover 
the entire range of interactions, from weakly interact- 
ing to strongly interacting. A few of the methods that 
have been applied to strongly-interacting systems include 
the coupled cluster method (GGM)[6], the method of 
correlated basis functions (GBF)[7], density functional 
theory [8], and quantum Monte-Garlo methods[5, 9, 10]. 

Dimensional perturbation theory (DPT) [11, 12] pro- 
vides a systematic approach to the study of correlation in 
quantum confined systems This method takes advantage 
of the high degree of symmetry possible among identical 
particles in higher dimensions to obtain an analytic de- 
scription of the confined quantum system - without mak- 
ing any assumptions about the number of particles or the 
strength of interparticle interactions. Because the per- 
turbation parameter is the inverse of the dimensionality 
of space (S = 1/D), DPT is equally applicable to weakly 
or strongly interacting systems. Another important ad- 
vantage of DPT is that low-order DPT calculations are 
essentially analytic in nature[13]. As a consequence the 



number of atoms enters into the calculations as a param- 
eter, and so, in principle, results for any A^ are obtainable 
from a single calculation [14]. Also, even the lowest-order 
result includes correlation, and so DPT may also be used 
to explore the transition between weakly-interacting sys- 
tems and those which are strongly interacting. These re- 
sults can be systematically improved by going to higher 
order [15]. 

In principle, excited states at low order in DPT are 
obtained from the same analytic calculation [14]. They 
only differ in the number of quanta in the different normal 
modes. Higher orders in the DPT expansion of excited 
states can also be calculated [15]. 

In past papers we have begun to apply dimensional 
perturbation theory to calculate the first-order ground 
state energy and normal mode frequencies of spherically- 
confined quantum systems[13, 14, 16-18]. In our last pa- 
per, we developed the formalism for deriving the lowest- 
order DPT wave function for an L = spherically con- 
fined quantum system (see Ref. [18] here also referred 
to as Paper I). In the present paper, we use the formal- 
ism developed in Paper I to calculate the lowest-order 
DPT density profile of the spherically confined quantum 
system. We apply these results to the specific case of a 
BEG for which the density profile is a directly observable 
manifestation of the quantized behavior of the confined 
quantum system. 

Section II provides a brief overview of the tools used. 
Section III discusses the L — 0, _D-dimensional, A^-body 
Schrodinger equation where L is the total angular mo- 
mentum of the system. Sections IV- VII provide a sum- 
mary of the relevant results from our previous work, par- 
ticularly from Paper I. Then in Section VIII the lowest- 
order density profile is derived from the lowest-order wave 
function of Paper I (summarized Sees. IV- VII), while Sec- 
tion IX discusses how it is possible to optimize our lowest- 
order result by minimizing the contribution from higher- 
order terms. Sections X and XI apply these results to the 
case of a BEG. Section X sets up the problem for a BEG, 
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assuming a hard-sphere interatomic potential at _D = 3 , 
while Section XI discusses the details of minimizing the 
contribution from higher-order terms. Section XII dis- 
cusses the results for a BEC. Section XIII presents our 
conclusions. 



II. TOOLBOX 

The tools used to describe large-iV correlated wave 
functions are carefully chosen to maximize the use of 
symmetry and minimize the dependence on numerical 
computation. We handle the massive number of inter- 
actions for large N (~ two-body interactions) by 
bringing together three theoretical methods. 

The first, DPT[11, 12], is chosen because its D ^ oo 
equation yields a maximally-symmetric configuration for 
N identical particles which allows an analytic solution. 
Higher orders yield insight into fundamental motions 
as well as a framework for successive approximations. 
The second method is the FG method of Wilson, De- 
cius, and Cross[19]. This seminal method has long been 
used in quantum chemistry to study vibrations of poly- 
atomic molecules. It directly relates the structure of 
the Schrodinger equation to the coordinate set which 
describes the normal modes of the system. The third 
method, the use of group theoretic techniques[19, 20], 
takes full advantage of the point-group symmetry in the 
D 00 limit. 



III. THE i3-DIMENSIONAL AT-BODY 
SCHRODINGER EQUATION 

The Schrodinger equation for an A^-body system of 
particles confined by a spherically symmetric potential 



with a two-body potential gij is 
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In a £)-dimensional Cartesian space 
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and Qij = Vint 
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are the single-particle Hamiltonian and the two-body in- 
teraction potential, respectively, Xii, is the f*'^ Cartesian 
component of the i*'' particle, and Vconf is the confining 
potential. 

A. The Effective S'-Wave Schrodinger Equation. 

It is desirable to transform from Cartesian to internal 
coordinates. A convenient internal coordinate system for 
confined spherically symmetric (L = 0) systems is 



{l<i< N), and 



7y = COs{0ij) = /riTj 



(4) 



{1 < i < j < N) , which are the D-dimensional scalar 
radii of the N particles from the center of the con- 
fining potential and the cosines "fij of the N{N — l)/2 
angles between the radial vectors. Under this coordinate 
change the effective S'-wave Schrodinger equation in these 
internal coordinates becomes 



E 



2m,; 



dn 



D-1 d 
n dn 



EE 



£)- 1 ^ d 



JV-l N 

Konf(ri) ^ + E E ^-t(^ii) 



(5) 



r 



where [n^f = [nf + [r^f - 2 7^^- . 

B. The Jacobian- Weighted Schrodinger Equation 

Dimensional perturbation theory utilizes a similarity 
transformation so that the kinetic energy operator is 

transformed into a sum of two types of terms, namely, 
derivative terms and a repulsive centrifugal-like term. 



The latter repulsive centrifugal-like term stabilizes the 
system against collapse in the large- D limit when at- 
tractive interparticle potentials are present. Low orders 
of the dimensional {1/D) expansion of the similarity- 
transformed Schrodinger equation are then exactly sol- 
uble for any value oi N . In the D ^ 00 limit the 
derivative terms drop out resulting in a static problem, 
while large-D corrections correspond to simple harmonic 
normal-mode oscillations about the infinite-dimensional 
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structure. (See Sections IV - VII below.) 

In Paper I the weight function was chosen to be the 
square root of the inverse of the Jacobian, J, where[21] 
J = (rir2...riv)^-^r(^-^-i)/2 and T is the Gramian 
determinant, the determinant of the matrix whose ele- 
ments are 7ij (see Appendix D of Ref. [13]), so that the 
similarity-transformed wave function $ and operators O 
are $ = J2 v[/ and O — O , respectively. Un- 
der this Jacobian transformation, a first derivative of an 
internal coordinate is the conjugate momentum to that 
coordinate. The matrix elements of coordinates and their 



derivatives between the lowest-order normal-mode func- 
tions, which are involved in the development of higher- 
order DPT expansions, are much easier to calculate since 
the weight function in the integrals is now unity. 

Carrying out the above Jacobian transformation of the 
Schrodinger equation of Eq. (1), we obtain [21]: 

{T + V)^ = E ^ (6) 

where 
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and 

JV N-l N 

v = j2 ^conf (n) + E E ) • (8) 

i=l i=l j=i-\-l 

Equation (7) for T is explicitly self-adjoint since the 
weight function, W , for the matrix elements is equal to 
unity. The similarity-transformed Hamiltonian for the 
energy eigenstate ^ is H = {T + V) . 

IV. INFINITE-I3 ANALYSIS: LEADING ORDER 
ENERGY 

Following Paper I, we begin the perturbation analysis 
by defining dimensionally scaled variables: 

fi = ri/K{D) , E = k{D)E , and H = k{D) H (9) 

where k{D) is a dimension-dependent scale factor, which 
regularizes the large-dimension limit. From Eq. (7) the 
kinetic energy T scales in the same way as l/r^ , so the 
scaled version of Eq. (6) becomes 

m= (^-^f + V,,,^^ = E^, (10) 

where barred quantities indicate that the variables are 
now in scaled units. The centrifugal-like term in T of 
Eq. (7) has quadratic D dependence so the scale fac- 
tor k{D) must also be quadratic in D , otherwise the 
_D — !■ 00 limit of the Hamiltonian would not be finite. 
The precise form of k{D) depends on the particular sys- 
tem and is chosen so that the result of the scaling is as 
simple as possible. In previous work[13] we have cho- 
sen k{D) = {D - 1){D - 2N - 1)/(4Z) for the ^-wave, 



A/'-electron atom; fll^o for the AT'-electron quantum dot 
where fi = {D—1){D — 2N— l)/4 and the dimensionally- 
scaled harmonic oscillator length and trap frequency re- 
spectively are ?ho \/^^ ^-^d = fi^w^^ ; and D^Sho 

for the BEC where Oho = \J^^^ and Who = -D^Who • In 
the 5 — > (Z) — > 00) limit, where 

5=1/D, (11) 

the factor of k{D) in the denominator of Eq. (10) sup- 
presses the derivative terms leaving behind a centrifugal- 
like term in an effective potential, 

Kff (r, 7; 5 = 0) =E y-^^ + ^--f r, 6 = 0) j 

iV-l N 

+ E E ^int(^,7;^ = o), 

i=l j=i+l 

(12) 

in which the particles become frozen. In the £> — > 00 

{6 0) limit, the excited states collapse onto the ground 
state at the minimum of Ve±f ■ 

We assume a totally symmetric, large-dimension con- 
figuration at which the effective potential is a minimum. 
The N particles are arranged on a hypersphere, each par- 
ticle with a radius, foo, from the center of the confining 
potential. Furthermore, the angle cosines between each 
pair of particles takes on the same value, 7^^, i.e. 

lim fi = foo i'i- < i < N) , 
hm 7y = 7^ (1 < i < i < iV) . 

D— >(X) 

This configuration for AT = 1, 2, 3, and 4 is comparable to 

the sequence of arrangements of hydrogen atoms in LiH, 
H2O, NH3 , and CII4 . In this analogy, the heterogeneous 
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atom represents the center of the confining potential [22]. 
(This symmetric high-dimensional structure is also not 
unlike the localized structure found in a hyperspherical 
treatment of the confined two-component normal Fermi 
gas in the N ^ oo Hmit[23].) 

In scaled units the D oc approximation for the en- 
ergy is simply the effective potential minimum, i.e. 



(14) 



In this D ^ oo approximation, the centrifugal-like term 
that appears in V^ft , which is nonzero even for the 
ground state, is a zero-point energy contribution satis- 
fying the minimum uncertainty principle [24]. 

Beyond-mean-field effects are already present in this 
approximation. This may be seen in the value of , 
the D 00 expectation value for the interparticle angle 
cosine (see Eqs. (13) and (17)). In the mean-field ap- 
proximation the expectation value for the interparticle 
angle cosine for the L = states considered in this paper 
is zero. The fact that 7,^ is not zero is an indication 
that beyond-mean-field effects are included even in the 
D 00 limit. 

This symmetric structure in which all N particles are 
equidistant and equiangular from every other particle can 
only exist in a higher-dimensional space and is impossible 
in a three-dimensional space unless TV < 4 (see above). 
That this high dimensional structure makes sense can 
be seen as follows. As we have noted above, foo and 
7;^ are the lowest-order DPT expectation values for fj 
and 7,j . This would indicate that the expectation values 
for the radii and interparticle angles for actual D = 3 
systems should have values corresponding to structures 
that can only exist in higher {D > 3) dimensional spaces. 
Accurate configuration interaction calculations for atoms 
in three dimensions do indeed have expectation values for 
the radii and inter-electron angles which define higher 
dimension structures [25]. 

This highly-symmetric, D ^ 00 structure imparts a 
point-group structure to the system which is isomorphic 
to the symmetric group of A'' identical objects[20], Sn, 
and allows for a largely analytic solution to the prob- 
lem, even though the number of degrees of freedom be- 
comes very large when N is large. The D ^ 00 ap- 
proximation may also be systematically improved by us- 
ing it as the starting point for a perturbation expansion 
(DPT) [15]. In this regard, the Sn symmetry greatly sim- 
plifies this task since the interaction terms individually 
have to transform as a scalar under the Sn point group. 



V. NORMAL-MODE ANALYSIS AND THE 1/D 
NEXT-ORDER ENERGY CORRECTION 

In the D ^ 00 limit, the particles arc frozen in 
a completely symmetric configuration (which is some- 
what analogous to the Lewis structure in chemical 
terminology[26]). This configuration determines the 
D —^ 00 energy, , but it is not a wave function as no 



nodal information is present. The large-dimension limit 
is a singular limit of the theory, so to obtain the lowest- 
order wave function we have to consider the next order 
in the perturbation expansion, which yields not only the 
(1/Z))°, i.e. leading-order wave function, but the order 
{l/Oy correction to the energy. The perturbation series 
has the form: 



00 

E(**)' 



00 



(15) 



3=0 



In practice Ej = y j odd. This next-order 1/D cor- 
rection in the energy, but lowest-order in the wave func- 
tion, can be viewed as involving small oscillations about 
the D ^ 00 structure (somewhat analogous to Lang- 
muir oscillations [2 7]). As we shall see below, these are 
obtained from a harmonic equation, and so we refer to 
them as the energy and wave function at harmonic order. 

To obtain this harmonic correction to the energy for 
large but finite values of £> , we expand about the mini- 
mum of the D ^ 00 effective potential. We first define a 
position vector, y , consisting of all A^(A^ -|- l)/2 internal 
coordinates: 



y 



7 



where 7 



and r 



\fNl 




(16) 



lN-2,N 
V IN-I.N I 



We then make the following substitutions for all radii and 
angle cosines: 

fi=foo+^'/M and 7u=7cx>+^'^% 



(17) 

In all practical situations l4ff (y; S) is a function of y and 

(5, so we may obtain a power series in b^l"^ of the effective 
potential about the D 00 symmetric minimum. 

Defining a displacement vector of the internal displace- 
ment coordinates [primed in Eqs. (17)] 



and 



where 7' 
1\ 




(18) 
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we find 
Veffiy'; S) 




where P = N{N + l)/2 is tlie number of internal coor- 
dinates. Tiie quantity, Vo, is just a constant term, in- 
dependent of y'^, while the first term of the 0((^^/^)^) 
term defines the elements of the Hessian matrix [2 8] F of 
Eq. (21) below. The derivative terms in the kinetic en- 
ergy are taken into account by a similar series expansion, 
beginning with a harmonic-order term that is bilinear in 
d/dy' , i.e. 



(20) 



where T is the derivative portion of the kinetic energy 
T see Eq. (7)). It follows from Eqs. (19) and (20) that 
G and F , both constant matrices, are defined in the 
harmonic-order Hamiltonian as follows: 



(21) 



Thus, obtaining the harmonic-order wave function, and 

energy correction is reduced to solving a harmonic equa- 
tion by finding the normal modes of the system. 



A. The FG Matrix Method for the Normal Modes 
and Frequencies 

We use the FG matrix method[19] to obtain the 
normal-mode vibrations and, thereby, the harmonic- 
order energy correction. A review of the FG matrix 
method is presented in Appendix A of Paper I, but the 
main results may be stated as follows. The fe**^ normal 
mode coordinate may be written as (Eq. (A9) of Paper I) 



(22) 



where the coefficient vector b satisfies the eigenvalue 
equation (Eq. (AlO) of Paper I) 



FGb^Xbb 



(23) 



with the resultant secular equation (Eq. (All) of Pa- 
per I) det(FG? — AI) = 0. The coefficient vector also 
satisfies the normalization condition (Eq. (A12) of Pa- 
per 1) b^Gb = 1. As can be seen from Eq. (A3) of 
Paper I the frequencies are given by Xb — iof . 

In an earlier papcr[13] we solve these equations for the 
frequencies. The number of roots A of the secular equa- 
tion - there are P = N{N + l)/2 roots - is potentially 
huge. However, due to the Sn symmetry of the prob- 
lem discussed in Sect. 5 and 6 of Paper I, there is a very 



significant simplification. The secular equation has only 
five distinct roots, A^ , where /i is a label which runs over 
0~ , 0+ , 1^ , 1+ , and 2 , regardless of the number of 
particles in the system (see Refs. [13] and Sect. VI B). 
Thus the energy through harmonic-order (see Eq. (25)) 
can be written in terms of the five distinct normal-mode 
vibrational frequencies which are related to the roots A^ 
of FG by 

>^^=^l- (24) 
The energy through harmonic order in 6 is then[13] 



E = Eoo+S 



E 

;*={0±,l=' 



n„=0 



2 } 



(25) 

where the are the vibrational quantum numbers of the 
normal modes of the same frequency o)^ (n^ counts the 
number of nodes in a given normal mode) . The quantity 
c^p.n^ is the occupancy of the manifold of normal modes 
with vibrational quantum number and normal mode 
frequency Cu^, i.e. it is the number of normal modes with 
the same frequency cS^ and the same number of quanta 
n^. The total occupancy of the normal modes with fre- 
quency (D^ is equal to the multiplicity of the root A^, i.e. 

= J2T=o '^M.n^ ) where is the multiphcity of the /x*'^ 
root. The multiplicities of the five roots are[13] 



do+ — 1 , do- = 1 , = — 1 , 
d^- =N-1, d2 = N{N - 3)/2 . 



(26) 



An analysis of the character of the normal modes re- 
veals that the 2 normal modes are phonon, i.e. compres- 
sional, modes; the 1^ modes show single-particle char- 
acter, and the 0^ normal modes in a harmonic con- 
fining potential describe center-of-mass and breathing 
motions. [29] 



VI. THE SYMMETRY OF THE F AND G 
MATRICES AND THE JACOBIAN- WEIGHTED 
WAVE FUNCTION 

A. Q matrices in terms of simple invariant 
submatrices 

Such a high degree of degeneracy of the frequencies 
in large- £) , A^-body quantum confinement problems in- 
dicates a very high degree of symmetry. The F , G , 
and FG matrices, which we generically denote by Q , 
are P x P matrices. The Sn symmetry of the Q matri- 
ces {F , G , and FG) allows us to write these matrices in 
terms of six simple submatrices which arc invariant under 
Sn . We first define the number of jij coordinates to be 
M = N{N - l)/2 , and let In he a,n N x N identity ma- 
trix, Im an M X M identity matrix, Jn an N x N matrix 
of ones and Jm an MxM matrix of ones. Further, we let 
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Rhe an N X M matrix such that R; 



- Sij + Sik , Jnm 
= Jmn ■ These 



be an X M matrix of ones, and i/j^^tm ^ *Jmn 
matrices are invariant under intercliangc of the particles, 
effected by the point group Sn- They also have a specific 
interpretation in the context of spectral graph theory (see 
Appendix B of Ref. [13]). 

We can then write the Q matrices as 



Q = 



Qi.j 


Qi.jk 


Qij,k 


Qij.kl 



(27) 



where the block Qfifi has dimension [N x N) , block 
QjHyi has dimension (NxM) , block Qy^/ has dimension 
{M X N) , and block Q^'^' has dimension (M x M) . As 
shown in Appendix B of Ref. [13], 



Qr'r' = {Qa — Qb)^N + Qb^N 1 
Qt='^' = {Qe — Qf)^ + Qf^NM 
Q^'fi = {Qc — Qd)R + <9dJ ATM 



(28) 
(29) 



Q^'t' = (Qs - 2Qfe + QJIm + (Qh - QOR^R 

+ g, Jm . (30) 

It is this structure that causes the remarkable reduction 
from a possible P = N{N + l)/2 distinct frequencies to 
just five distinct frequencies for L — systems[30]. 

In particular, letting Q be FG , the matrix to be diag- 
onalized, or G , the matrix required for the normalization 
condition, Eq. (27) becomes 



FG 



or G = 



oIn + bJN 
cR^ + dJ 



MN 



eR + f^NM 
glu + hR^H + 13 



M 



(31) 



g'lM + /I'R^R 



where the coefficients a , b , c , d , e , f , g , h , and t have a 
simple relation to the elements of the F and G matrices 
while d' , g' , and h' have a simple relation to the elements 
of the G matrix (see Paper I). These coefficients also 
depend on the choice of k{D) (See Ref. [13]). 



B. Symmetry Coordinates, Normal Coordinates, 
and the Jacobian-Weighted Wave function 

As discussed in Paper I the FG matrix is a N{N + 

l)/2 X N{N + l)/2 dimensional matrix (there being 
N{N + l)/2 internal coordinates), and so the secular 
equation could have up to A''(A'' + l)/2 distinct frequen- 
cies. However, as noted above, there are only five dis- 
tinct frequencies. The Sn symmetry is responsible for 
the remarkable reduction from N{N + l)/2 possible dis- 
tinct frequencies to five actual distinct frequencies. As 
we shall also see, the Sn symmetry greatly simplifies the 
determination of the normal coordinates and hence the 
solution of the large-D problem. 

The Q matrices, and in particular the FG matrix, 
are invariants under Sn , so they do not connect sub- 
spaces belonging to different irreducible representations 



of ^^[Sl]. Thus from Eqs. (22) and (23) the normal co- 
ordinates must transform under irreducible representa- 
tions of Sn ■ Since the normal coordinates will be linear 
combinations of the elements of the internal coordinate 
displacement vectors f' and 7' , we first look at the Sn 
transformation properties of the internal coordinates. 

The internal coordinate displacement vectors f' and 7' 
of Eqs. (17) are basis functions which transform under 
matrix representations of Sn , and each span the corre- 
sponding carrier spaces, however these representations of 
Sn are not irreducible representations of Sn ■ 

In Sec. 6.1 of Paper I we have shown that the re- 
ducible representation under which r' transforms is re- 
ducible to one 1-dimensional irreducible representation 
labelled by the partition [N] (the partition denotes a cor- 
responding Young diagram ( = Young pattern = Young 
shape)) of an irreducible representation (see Appendix C 
of Paper I) and one {N — l)-dimensional irreducible rep- 
resentation labelled by the partition [N — 1, 1] . We 
will also show that the reducible representation under 
which 7' transforms is reducible to one 1-dimensional ir- 
reducible representation labelled by the partition [N] , 
one (A'^ — l)-dimcnsional irreducible representation la- 
belled by the partition [TV - 1, 1] , and one N{N - 3)/2- 
dimensional irreducible representation labelled by the 
partition [A'' — 2, 2] . Thus if da is the dimensionality 
of the irreducible representation of Sn denoted by the 
partition a then d[N] = 1 , '^[w-i. 1] = — 1 , and 
d[N-2, 2] = N{N—3)/2 . We note that: d[N]+d[N-i, i] = 
AT , giving correctly the dimension of the f' vector and 
that: d[N] + cf[w-i, 1] + rf[7V-2, 2] = N{N - l)/2, giving 
correctly the dimension of the 7' vector. 

Since the normal modes transform under irreducible 
representations of Sn and are composed of linear com- 
binations of the elements of the internal coordinate 
displacement vectors f' and 7' , there will be two 1- 
dimensional irreducible representations labelled by the 
partition [A''] , two (A'' — l)-dimensional irreducible repre- 
sentations labelled by the partition [A^ — 1, 1] , and one 
entirely angular A^(A^ — 3)/2-dimensional irreducible rep- 
resentation labelled by the partition [A^ — 2, 2] . Thus if 
we look at Eq. (26) we see that the 0* normal modes 
transform under two [A^] irreducible representations, the 
1=*= normal modes transform under two [N — 1, 1] irre- 
ducible representations, while the 2 normal modes trans- 
form under the [A^ — 2, 2] irreducible representation. 

The wave function for the harmonic-order Hamilto- 
nian of Eq. (21) is thus the product oi P = N{N + l)/2 
harmonic-oscillator wave functions: 



*o(y') = n n 

/:i={0±,l±,2} £=1 



<l>nJV^ [q'\) , (32) 



where (pn^^^^/uF^ W^]^) is a one-dimensional harmonic- 
oscillator wave function of frequency u),,, and n„,e is the 



which counts 



oscillator quantum number, < n^^ < 00 

the number of quanta in each normal mode. The quantity 
jj, labels the manifold of normal modes with the same 
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frequency while = 1 , A'' — 1 or N{N — 3)/2 for 
/U = 0^ , l"*" or 2 respectively. 



VII. CALCULATING THE NORMAL MODE 
COORDINATES. 

In Paper I, we extended previous work to the calcula- 
tion of the normal coordinates. We summarize this work 
below in two-steps: 

a) . In Paper I we have determined two sets of 

linear combinations of the elements of coor- 
dinate vc;c;tor f' which transform under par- 
ticular orthogonal [N] and [A^ — 1, 1] irre- 
ducible representations of Sn ■ These are 
the symmetry coordinates for the f' sector 
of the problem [19]. Using these two sets of 
coordinates we then determined two sets of 
linear combinations of the elements of coor- 
dinate vector 7' which transform under ex- 
actly the same orthogonal [N] and [A'' — 1, 1] 
irreducible representations of Sn as the co- 
ordinate sets in the r' sector was deter- 
mined. Then another set of linear combi- 
nations of the elements of coordinate vector 
7' , which transforms under a particular or- 
thogonal [A'' — 2, 2] irreducible representation 
of ^jv , was determined. These are then the 
symmetry coordinates for the 7' sector of the 
problcm[19]. Furthermore, we chose one of 
the symmetry coordinates to have the sim- 
plest functional form possible under the re- 
quirement that it transforms irreducibly un- 
der Sn ■ The succeeding symmetry coordi- 
nate was then chosen to have the next sim- 
plest functional form possible under the re- 
quirement that it transforms irreducibly un- 
der Sn , and so on. In this way the complex- 
ity of the functional form of the symmetry 
coordinates has been kept to a minimum and 
only builds up slowly as more symmetry co- 
ordinates of a given species were considered. 

b) . The FG matrix is expressed in the f ' , 7' 

basis. However, if we change the basis in 

which the FG matrix is expressed to the sym- 
metry coordinates an enormous simplification 
occurs. The N{N +l)/2 x N{N +l)/2 eigen- 
value equation of Eq. (23) is reduced to one 
2x2 eigenvalue equation for the [A^] sector, 
A^ — 1 identical 2x2 eigenvalue equations for 
the [A^-1, 1] sector and A'(A'-3)/2 identical 
1x1 eigenvalue equations for the [A'' — 2, 2] 
sector. In the case of the 2x2 eigenvalue 
equations for the [A^ - 1, 1] and [AT - 2, 2] 



sectors, the 2x2 structure allows for the mix- 
ing in the normal coordinates of the symme- 
try coordinates in the f' and 7' sectors. The 
1x1 structure of the eigenvalue equations 
in the [A' — 2, 2] sector reflects the fact that 
there are no [A^ — 2, 2] symmetry coordinates 
in the f' sector for the [A'' — 2, 2] symmetry 
coordinates in the 7' sector to couple with. 
The [N — 2, 2] normal modes are entirely an- 
gular. 



A. Transformation to Symmetry Coordinates 



Let the symmetry coordinate vector, S , be 



aim 



gjN-l, 1] 



S[N] 

= I 5'[A'-1, 1] 
S[N-2, 2] 



(33) 



In S, symmetry coordinates of the same species are 
grouped together. Prom Paper I we have 



[4^ 



N{N - 1) 



JV 



(34) 



j=2 i<j 



1. 111. _ 



where 1 < i < N ■ 



ir. 



1 



(35) 



We see from Eq. (35) the first symmetry coordinate is 
proportional to — . This involves only the first two 
particles in the simplest motion possible under the re- 
quirement that the symmetry coordinate transforms ir- 
reducibly under the [A^ — 1, 1] representation of Sn ■ We 
again advert that adding another particle to the system 
does not cause widespread disruption to the symmetry 
coordinates. The symmetry coordinates, and the motions 
they describe, remain the same except for an additional 
symmetry coordinate which describes a motion involving 
all of the particles. 

Again according to item a), above, S^, ^' ^' should 
transform under exactly the same non-orthogonal irre- 
ducible [A^ — 1, 1] representation of Sn as S^^~'^' . 
From Paper I 



where 1 < i < N — 1 . Now there is only one sector, Paper I 
the 7' sector, belonging to the [N — 2, 2] species. Prom 



[N-2, 2h 



Vi(i + l)(i-3)(i-2) 



j-l [j'-l,i]„i„ i-1 i i 

E E + E E ^^i' - (-^^ - 3) E ^^i- 

j'=2 fc=l fe=l j'=fc+l fc=l 

i j-1 

-*E^M»+i)"^ E ^(i+i) J' + - 3)7(i+i) J 

fe=l j'=i+2 

I 



(37) 



where 1 < i < j - 2 and 4: < j < N . 

This fulfills item a), above: the determination of the 
symmetry coordinates. 



B. Transformation to Normal-Mode Coordinates. 

In Paper I we apply the FG method using these sym- 
metry coordinates to determine the frequencies and nor- 
mal modes of the system: 



A„ — 



(38) 



a„ = 



[CTa \f',f' 



'^a.^\r',y [fa JT'.r' ) 



for the a = [N] and [N — 1, 1] sectors, where 

(39) 

while A[;v-2, 2] = '^[n'—^, 2] • '^^^ related to 

the 5, , & , c , d , e , / , g , h , and l of the i^G matrix 
of Eq. (31). Note that the cr^*^ for the a = [N] and 
[A'' — 1, 1] sectors are 2x2 matrices while o'^_2^ 2] 
a one-component quantity. 

The normal coordinates are given by 

q'l = cl (cos ei [S^,]^ + sin [S^,]^) (40) 

for the a = [N] and [A'' — 1, 1] sectors, and 

^,[JV-2, 2] ^ ^[jV-2, 2] _5^W-2, 2] ^^^^ 



The r'-7' mixing angle, 9^. , is given by 



tan6'i = 



(AJ-[T^^].',.') 



(42) 



while the normalization constants c are given by 
1 



c± — 



and 



sm( 



uub iy_|_ 

sin^^ 



(43) 



4JV-2, 2] _ 



\J'^fN-2, 2] 



The cr^ are related to the a' , g' , and /i' of the G matrix 
of Eq. (31). 

This fulfills item b). above; the determination of the 
frequencies and normal coordinates. Equations (40) and 
(41) represent the final results of Paper I, the normal 
coordinates. 



VIII. DERIVING THE DENSITY PROFILE OF 
A BEC FROM THE WAVE FUNCTION 



As a first application of this general theory from Pa- 
per I, we calculate the density profile for a confined quan- 
tum system. The harmonic-order DPT wave function, 
g$o(y') ) for the ground state is given by Eq. (32) with 
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all of the n^^ set equal to zero, i.e. 



My') = n n ' (44) 

,i={0±, 1±, 2} 1=1 

where 

HV^. = (^) " exp . (45) 

Now the ground-state wave function of a BEC is com- 
pletely symmetric under interchange of any of the par- 
ticles and so we require g^o{y') of Eq. (44) to be 
completely symmetric under interchange of any of the 
particles. This follows from the fact that the five 

nj=i Ma/^ W^k) ' where ^ = {0="=, 2} , are each 
completely symmetric under interchange of any of the 
particles. This may be seen from the fact that 




E } I . (46) 



Now according to items a), and b). of Sec. VII, particle 
interchange is effected by orthogonal transformations of 
the [q'^]^ , where 1 < ^ < , and any orthogonal trans- 
formation leaves X^^^i invariant. Together with 

Eq. (46), this fact means that the Yl'^'Li 4>a{\^^ W^]i) 
arc completely symmetric under interchange of any of the 
particles as advertised. 

Defining S{D) to be the total £)-dimensional solid 
angle [32], 



S{D) 



2 TT^ 

f(¥)' 



(47) 



where we note that 5(1) = 2, 5(2) = 27r, 5(3) = An, 
5(4) = 27r^ , • • • , the harmonic-order Jacobian-weighted 

ground-state density profile, Noir) , is 

S{D)No{r) = S{D)r^"-^^po{r) 

= E / • • • / ^f(^ - ^0 [«^o(y')]' n n ' 

i=l-^-°° -^-oc ^=0±,1±,2€=1 

(48) 

where po{r) is the unweighted harmonic-order ground- 
state density profile and (5/(r — Tj) is the Dirac delta 
function and is differentiated from the inverse dimension, 
S, by the subscript /. Since g$o(y') is invariant under 
particle interchange then 

S{D)No{r) 

.../ Sf{r-rj,)[My')r n n^[q'i«- 

-<X J — CO .._„+ „ c_i 

(49) 



/:i=0±,l±,2 5=1 



Upon using Eq. (44) and the fact that rjv only appears in 

q'°^ , q'°~ , [q'^^]d^+ and [q'^'h (see Eqs. (9), (17), 
(34)-(37), (40), (41)) and Eq. (45)) we obtain 



S{D)No{r) 

pOO pOD pOD pOO 

N ^ 

-OO J — OO J — OO J — OO 

/n=0±,l± 



/OO pOO pOO pOO 
/ / / Sf{r-rN) 
-OO J — OO J — OO J — OO 



Ny/CUQ+ UJq^ UJi+ LUi- 



/OO pCO pC<D pOO 
/ / / 6f{r-rN) 
-OO J — OO J — OO J — OO 



exp I - ^ co.iqni I n dW%c 



. i/=0±,l± 



C=o±,i± 



(50) 

The delta function, 5f{r — tn) , is a function of rjv 
while the integral is over the normal coordinates q' , 

9'° I ]n-i , and [q'^ ]n-i ■ Thus we need a 

change of variables to perform the integral. We change 
the variables of the integral to f'f 



N ' 



, 5[^J, and 



[N-l, 1] 



J(W-l) 



, where 



N-l 



E^^ 



(51) 



SL, is given by Eq. (34) and [S^' ^' ^'](Ar_i) is given by 
Eq. (36). Thus from Eqs. (35), (40), (41), and (51) we 
obtain 



( q'^^ \ 



\N-1 



= T a' , 



(52) 



V W \n-i J 



where 



' N 



5^1 



I tqlN-l, l]l / 

\ [^^' \{N-1) I 



(53) 



/ Tn Tn T13 \ 

Tlx T21 T23 

-{N-l)Tsi T31 T34 

y _(7V-i)T4i T41 T44 / 



and 
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T31 = 



N 

[N-l, 1] 



Ti3 = cflsin0fJ 



To, = 



[N] alN] 

c_ ' COS 0^ ' 



\/iV 



, -123 — C_ falll(7_ 



^[Af-l, 1] 



V^(iv - 1) 



-t 34 — C_,_ 



. /,[JV-1, 1] 

sin 6*^ ' 



(54) 



T44 = cL^-^' ^Isin^L^-^' ^1. 

I 



The Jacobian, Jt , of the transformation is thus 



Jt = det T 



JN] JN] [N-l, 1] [JV-1, 1] 
O j C- -|- 



X sin(^fl -^L^]) sin(^f-^' 1] -eL'^-i- 11). (55) 



Defining 



Wo+ 

Wo- 

^1+ 

wi- 



(56) 



then 

S{D)No{r) 

N Jt^/ojo+ ojq- a>i- 



00 /'OO /'CO /'OO 



00 00 -/—CO 00 



(5/(r-rjv) 



where d^b' = dr'sdS^^"^ d[S^^~^' \n-i) ■ Writing 



T^nT 



1\0 J\' 

K K 



(58) 



where 



(57) 



(59) 



/ iOo+Tfi + iOo-Tl -{N- l)(c^i+T|i + Wi-Tli) 
-K" = W0+T11T13 + a;o-T'2iT23 

\ -{N - l)(a;i+r3iT34 + ujx-TiiTa) 



/Cii 
/Ci3 = Ksi 
/C23 = /C32 



= Wo+Tfi + Wo- ^21 + '^1+^31 + ^i-T^i , K\2 = K.21 = a;o+TiiTi3 + a;o-^2i75 



23 : 



'^l+T'3lT34 + UJi-TiiT^i , K22 = + Wo- ^23 ) 

0, K33=u)i+Tl+u)i-T^ 



(60) 



(61) 



144 • 



Using Eqs. (53) and (58) in Eq. (57), we obtain 

S{D)N„{r) 

= ^5 / / / / ^/(^-^^) 

X exp (-iiTor^ - If'^K'^b' - b'^Kb') df'^ d^b' , 

(62) 



where 



b' = 



\ 



(63) 



Since df{r-rN) = VS k{D) Sf (d 2 (^-^ - f^'j - 
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then 



S{D)No{r) = 



N JtVS k{D) ^/LOo+ LOq- LOi+ LOi- 



/oo |.oo |.oo |.oo / iff 

.00 J -00 J -00 7-00 v "'^'^ v ^ 

26- i I -4^ -rS] K^b' - b'^Kb'] d?b' 



(64) 



Using 



/OO nOC 
-00 J — 



exp (-b'^Ab' - 2B^6') d"b' 



00 */ —00 



7r2 



Vdet A 



exp(B^A-iB) 



(65) 



Eq. (64) yields 



X exp i —6 ^ ( — ^ 



V det /C 

2 



(66) 



Since 



— — — — t2 
'^0+ ''^O- ^1^1+ Wi- 

detX; 



= (ifo-ii:'/C-^K), (67) 



then 



S{D)No{r)^N- 



X exp I ^ ( 



(68) 



The harmonic-order DPT Jacobian-weighted den- 
sity profile, No{r) , is correctly normalized to N, 
since upon using Eq. (65) (with n = 1 , A = 
{Ko - K'^K-^K)/{6k^{D)) and B = -f^{KQ - 
K'^K-^K) {6k{D))) in Eq. (68) we find 



/ / No{r) dr dn = N . 



(69) 



Notice that the harmonic-order DPT density profile is a 
normalized gaussian centered around r = K{D)fao , the 
2? — > 00 configuration radius (see Eqs. (9) and (13)). 



IX. OPTIMIZING THE HARMONIC-ORDER 

RESULT 

Let's look in more detail at the harmonic-order wave 
function. In the notation of Eqs. (16) and (18), Eq. (17) 



can be written as 



where 



y = yoo 



Tjfc Toe 



+ 5^'^y', (70) 



V 1 < i < and (71) 
l<i<k<N. 



Inserting Eq. (70) into Eq. (22) one obtains q — Qo 



S^/'^q' , where Qo 
obtains 



b Voo ■ Using these in Eq. (32) one 



*o(y') = n n <i^n,,{{Du^,Y'' mi - g^)) . 

;i={0±, 1±, 2} «=1 

(72) 

Equation (72) represents oscillations about the Z) — > 00 
configuration qr^ with frequencies {Du)^} . For a macro- 
scopic quantum confined system at i5 = 1/3 {D = 3), 
Dojfj, = Sojfj, is sufficiently small that this harmonic- 
order, Jacobian-weighted wave function has a macro- 
scopic extension. However, when D becomes large {6 be- 
comes small) the frequencies {Duj^^} becomes very large 
and so, according to Eq. (72), the harmonic-order wave 
function becomes strongly localized about [q^]^ = qf^ 
(i.e., it features short wavelengths). 

This dimensional behavior has interesting, and very 
useful consequences. As is well known, the energy and 
density profile of a small-a Bose-Einstein condensate at 
D = 3 depends only on the scattering length of the in- 
teratomic potential, and not the detailed shape of the 
potential. This is due to the long wavelength nature 
of EEC's: for small to moderate scattering lengths, the 
atomic wavelengths are not short enough to "resolve" the 
short-range detail of the potential. However, for large D 
the atomic wavelengths become very short, since accord- 
ing to Eqs. (9) and (10), the scaled, Jacobian-transformed 
Hamiltonian displays an effective mass term proportional 
to . Thus as we have noted above, and unlike at 
D = 3 , the wave length of the wave function for the 
large-D system becomes smaller and becomes sensitive 
to the details of the potential. 

This feature is actually advantageous. A perturbation 
theory in some parameter which at low orders displays 
an insensitivity to the precise shape of the interatomic 
potential could not be optimized to yield results at low 
order that would be close to the actual results. For ex- 
ample, one could not reasonably expect the energy and 
density profile at low orders to be both insensitive to the 
precise shape of the interatomic potential for small fixed 
scattering length and, at the same time, to differ only a 
small amount from the actual D = 3 condensate energy 
and density profile. The energy and density profile at 
low orders would almost certainly be different from the 
actual D = 3 condensate. 

The large- Z? sensitivity to the details of the interatomic 
potential in the present method enables us to optimize 
our dimensional continuation of the interparticle poten- 
tial so that higher-order contributions of this theory will 
be small. 
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X. APPLICATION: THE DENSITY PROFILE 
FOR AN AT-ATOM BOSE-EINSTEIN 
CONDENSATE 

We assume a T = OK condensate which is confined by 
an isotropic, harmonic trap with frequency uiho- 



(73) 



A realistic two-body atom-atom interaction potential 
would lead to a solid-like ground state, and so typically 
this two-body potential is replaced with a pseudopoten- 
tial which does not support bound states. We follow most 
other work and replace the actual atom-atom potential 
by a hard sphere of radius a: 



DO , 

0, 



rij > a . 



(74) 



where a is the s-wave scattering length of the conden- 
sate atoms. We dimensionally continue the hard-sphere 
potential so that it is differentiable away from £) = 3 , al- 
lowing us to perform the dimensional perturbation anal- 
ysis (see Ref. [13, 14] as well as a later discussion in this 
paper). Thus, we take the interaction to be 



Vo 



1-3/D 



1 — tanh 



1-3/D 



{rij - a) 



(75) 

where D is the Cartesian dimensionality of space. This 
interaction becomes a hard sphere of radius a in the phys- 
ical, D = 3 , limit. The two constants, Vo and Cq (which 
determine the height and the slope of the potential), 
arc parameters that allow us to fine-tune the largc-_D 
shape of the potential and optimize our results through 
harmonic order by minimizing the contribution of the 
higher-order terms beyond harmonic (see Section XI). 
Although we have chosen the simplest possibility for the 
interatomic potential, two parameters, we can have any 
number of parameters providing for a more general and 
flexible potential[14]. The functional form of the poten- 
tial at 7^ 3 is not unique. Other functional forms could 
be chosen with equal success as long as the form is differ- 
entiable and reduces to a hard-sphere potential at £> = 3 . 
We simply choose a form that allows a gradual softening 
of the hard wall. 

Wc need to regularize the large-D limit of the 
Jacobian-weighted Hamiltonian (J^/^iJ J~^/^). We do 
this by converting the variables to dimensionally-scaled 
harmonic-oscillator units (bars): 



= Co = V2D^ahoCo , 

where 



V2D^aho' 

(76) 



dho 



mtoho 



and ujho = D^Who (77) 



are the dimcnsionally-scalcd harmonic-oscillator length 
and dimcnsionally-scalcd trap frequency, respectively. In 
this case wc have chosen k{D) = D^a^a , while the 
dimensionally-scaled harmonic-oscillator units of energy, 
length and time are fwoho , auo , and 1 /uho , respectively. 
All barred constants (a , aho , i^ho , Vo , and Co) are held 
fixed as D varies. For example, as D varies a is held fixed 
at a value by requiring that it give the physical unsealed 
scattering length at Z) = 3 . 

In dimensionally scaled units the total interaction term 
reads 



V = Ko.f + V,, 



^ 1 



2 ' ^ 1 



Vo 



N 



35 



N 



1 — tanh 



1-3^ V\/2 



— a 



.(78) 



The infinite-Z) {6 
dimensionally-scaled 
Eqs. (76) and (77) is 
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0) effective potential in 
harmonic-oscillator units of 




=1 



V2 



■ (79) 



One can see from the double-sum term in V^ff that the 
large-D interatomic potential has become a soft sphere 
of radius a and height 2Vo ■ The slope of the soft wall 
is determined by Cq . The development of DPT using 
the three-parameter extension of Eq. (75) is discussed at 
length in Ref. [13], while a four parameter extension of 
Eq. (75) is discussed in Ref. [14]. 

As noted above, the two parameters in Eq. (79) are 
chosen with the goal of minimizing the contribution of the 
higher-order beyond-harmonic terms, and so optimizing 
the harmonic-order DPT density profile (of Eq. (68)) and 
energy perturbation series through harmonic order in 6 
(Eq. (25)). 



XI. OPTIMIZATION OF THE INTERPARTICLE 
PARAMETERS OF A QUANTUM CONFINED 
SYSTEM 

We consider a BEG in a spherical trap where LUho = 
2-K X 77.87Hz and for which a = 1 OOP a.u. or OMSSaho 
in oscillator units {aho = y^fi/mtOho) which is approx- 
imately equal to ten times the natural ^''Rb value. We 
choose this value for the scattering length since the actual 
density profiles and energies show a noticeable difference 
from the mean-field, GP equation result at low N , and 
yet the modified Gross-Pitaevskii (MOP) equation[33] re- 
mains valid for comparison to the DPT result. 

The potential is optimized by fitting the energies 
through harmonic order to benchmark DMC energies[10] 
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at low atom number (A'' < 50) and DMC density pro- 
files for = 3 and N = 10 (DMC density profiles for 
larger TV have not been obtained[10]). Since in our DPT 
analysis the number of atoms iV is a parameter, we can 
readily extrapolate to larger N without large amounts of 
calculation. 

A least-squares fit is used to optimize the parameters 
of the dimensionally-continued interatomic potential. We 
fit to six accurate low-iV DMC cncrgics[10] and, as noted 
above, the position and height of the peak of the DMC 
density profile for A'' = 3 and N = 10 . The DPT values 
for the position and height of the peak of the density 
profile are just 



dDPT) 
peak 



{N;Vo,Co) = fo 



(80) 



S{D) V TT 



(81) 



D=3 



Thus we minimize the following quantity with respect 
to the two parameters Vo and cq: 



6 

i=l 



E, 



{DMC) 



E 

JV={3,10} 



+ 



{N) 



+ 



-£(^^^)(iv,;v;,c„)) 

(Ar(Sf^)(Ar) - [N,)^^:J\n-%,-Co)) 







-{DMC) 
peak 



(82) 

where e[^'^'~^^ is the dimcnsionally-scalcd DMC energy 
for a condensate with atom number Ni , while r^^^^^*^' {N) 

and Nmax'^\N) are the dimensionally-scaled position 
and height of the peak of the DMC density profile. The 
quantity ^(^^"^^ (iVj; V;, c^) is the DPT energy approxi- 
mation through harmonic order given by Eqs. (25) with 
interatomic potential parameters Vo and Cq , evaluated at 
£1 = 3. Thus we have a two-parameter fit for Co and Vo ■ 
This results in values Co = 4.7315 and Vo = 630.9573 at 
an i?2 of 14.6015. 

We remind the reader that the potential at D = 3 
remains a simple hard-sphere potential at r = a. 



XII. RESULTS AND DISCUSSION 

In Fig. 1 we plot the Jacobian-weighted density profile 
divided by the number of atoms N ai N = 100 for a 
scattering length a/aho ~ 0.0433, or roughly ten times 
the natural scattering length of ^''Rb when u>ho = 27r x 
77.87Hz . One himdrcd atoms is a ten-fold increase over 
the largest N density profile initially used to determine 
the interatomic potential parameters Co and Vo ■ As we 
see, the DPT density profile lies close to the MGP result, 
particularly in the 0.5 < r/uho < 1 range, although the 
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FIG. 1: The harmonic-order number density per atom versus 
radial distance of a spherically confined BEG of 100 *'^Rb 
atoms with a = 1 000 a.u. and ujho = 2tt x 77.87 Hz. The 
solid line is the analytic DPT density and the dashed line is 
the MGP density. Both curves are weighted by the Jacobian. 
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FIG. 2: The MGP number density per atom versus radial 
distance of a spherically confined BEG of 1,000 *^Rb atoms 
with a = 1 000 a.u. and coho — 2-kx 77.87 Hz. The dashed line 
is the MGP density weighted by the Jacobian. The vertical 
line is inserted to emphasize the asymmetry of the MGP curve 
at higher N . 



peak of the DPT profile is a little higher. Since the area 
under the curve is normalized to one, the extra height of 
the DPT peak comes at the expense of density at larger 
values of r/uho which lies below both the GP and MGP 
result. 

At larger values of TV the density profile develops an 
asymmetric aspect. This is illustrated in Fig. (2) where 
we plot the MGP Jacobian-weighted density profile at the 

same scattering length for 1, 000 . The lowest-order 
DPT density profile of Eq. (68) is a gaussian and symmet- 
ric about r^l^I'J'^ [N; Vo, Co) , and so is unable to capture 



this emerging asymmetry as A^ increases. Also if we fix 
A' and increase the scattering length yet further, the den- 
sity profile becomes increasingly asymmetric. Thus while 
the harmonic-order DPT density profile is reasonably ac- 
curate for small a and for smaller A^ at intermediate a, 
incorporating beyond- mean-field effects for more strongly 
interacting systems (larger TV or a) means we need to go 
to next order beyond harmonic in the DPT perturbation 
expansion to achieve an asymmetric beyond-mean-field 
Jacobian-weighted density profile. 

A next-order calculation for more strongly-interacting 
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systems may be outlined as follows. The higher- 
order beyond-harmonic interaction terms have to be re- 
expressed in terms of the normal coordinates of the large- 
D system. In this regard the Sn symmetry greatly sim- 
plifies this task since the interaction terms individually 
have to transform under the scalar [N] irreducible rep- 
resentation of Sn ■ In particular, there are only a finite 
number of Clebsch- Gordon coefficients of Sn coupling 
a finite number of normal coordinates, which transform 
under the three above discussed irreducible representa- 
tions of Sn , together to form a scalar [N] irrep. For ex- 
ample, there are only eight Clebsch-Gordon coefficients 
of Sn which couple three normal coordinates transform- 
ing under the [N] , [N - 1, 1] , or [TV - 2, 2] irreducible 
representations together to form a scalar [N] irrep. Ac- 
tually, it is this limited number of Clebsch-Gordon co- 
efficients which has allowed for the essentially analytic 
solution discussed in this paper of the harmonic-order 
DPT approximation and it greatly simplifies the higher- 
order beyond-harmonic interaction terms and their neces- 
sary transformation to normal coordinates. The required 
Clebsch-Gordon coefficients may be calculated analyti- 
cally for arbitrary N . 

XIII. CONCLUSIONS. 

In this paper wo have calculated an iV-body, analytic, 
lowest-order DPT Jacobian-weighted density profile for a 
quantum confined system with a general two-body inter- 
action from the previously derived analytic lowest-order 
DPT wave function[18]. The density profile is directly 
observable for macroscopic quantum confined systems, 
such as a BEC. 

The theory applied in this paper for L = states of 
spherically confined systems is applicable to systems with 
attractive or repulsive intcrparticlc interactions and is 
also applicable to both weakly (mean field) and strongly 
(beyond-mean-field) correlated systems. As we have 
noted above, the fact that 7,^ is not zero is an indication 



that beyond-mean-field effects are included in this result 
even in the D 00 limit. This theory is readily gener- 
alizable to systems for which the confining potential has 
cylindrical symmetry, 

We have applied the general formalism for the 
harmonic-order DPT Jacobian-weighted density profile 
developed in this paper, to the example of a spherically 
confined BEC. The higher-order DPT terms are mini- 
mized by analytically continuing the interatomic poten- 
tial in D away from a hard sphere of radius a at £> = 3 
so that we fit closely to accurate low-N DMC energies 
and Jacobian-weighted density profiles (N = 3 and 10 
for the density profiles). We have chosen a scattering 
length a/ttho = 0.0433, or roughly ten times the natu- 
ral scattering length of *^'^Rb when ojho = 27r x 77.87Hz . 
At this value of a the bcyond-mcan-field MGP equation 
remains valid for comparison. Exploiting the fact that 
our result is largely analytic, so that the particle num- 
ber, iV , is a parameter in our theory, we have tested the 
Jacobian-weighted density profile at = 100 and find 
the DPT result to lie close to the MGP result. 

However, for more strongly interacting systems the 
Jacobian-weighted density profile develops an asymmetry 
about the peak that our harmonic-order DPT Jacobian- 
weighted density profile does not mimic. Thus higher- 
order calculations arc required for larger N and large 
a . A detailed program for calculating higher-order DPT 
corrections to A''-body systems has been laid out in the 
paper by Dunn et. al.[15] and has been applied to high 
order for small- A^ systems [34]. As discussed at the end 
of Section XII of this paper, for large- A^ systems the Sn 
point-group symmetry is at the heart of, and greatly sim- 
plifies, the calculation of these higher-order terms. 
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